{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 量子相位处理\n",
    "\n",
    "*Copyright (c) 2022 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概览"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**量子相位处理**（quantum phase processing，简称 **QPP**）是一个由百度量子团队 [[1]](https://arxiv.org/abs/2209.14278) 提出的量子算法框架。该框架可以转换和处理酉算子的本征相位，并高效且精准地进行酉算子的本征变换或信息提取。该框架的推出来源于一种名为**三角量子信号处理**（trigonometric quantum signal processing, 简称三角 QSP）技术的改进 [[2]](https://arxiv.org/abs/2205.07848)。该技术可以使用单个比特去模拟输入数据的任意三角变换。对应地，QPP 则在更高维度上继承了三角 QSP 的能力。QPP 通过控制酉算子的方式，获取了酉算子的本征信息，从而在其本征空间下完成了对应本征相位的任意三角变换。\n",
    "\n",
    "基于论文 [[1]](https://arxiv.org/abs/2209.14278) 和块编码（block encoding）技术，本教程会向大家介绍如何使用 QPP 去解决量子相位估计问题、哈密顿量模拟问题和熵估计问题。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下是必要的 libraries 和 packages。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import numpy as np\n",
    "from scipy.linalg import expm\n",
    "import paddle\n",
    "\n",
    "# 量桨的通用函数\n",
    "import paddle_quantum as pq\n",
    "from paddle_quantum.ansatz import Circuit\n",
    "from paddle_quantum.hamiltonian import Hamiltonian\n",
    "from paddle_quantum.linalg import abs_norm, hermitian_random, unitary_random, block_enc_herm, herm_transform\n",
    "from paddle_quantum.qinfo import trace_distance, partial_trace_discontiguous\n",
    "from paddle_quantum.state import is_density_matrix, random_state, zero_state, State\n",
    "\n",
    "# 量桨的 QPP 模块函数\n",
    "from paddle_quantum.qpp import Q_generation, hamiltonian_laurent, qpp_angle_approximator, qpp_cir, qps, qubitize, purification_block_enc, simulation_cir\n",
    "\n",
    "# 设置计算精度和后端\n",
    "pq.set_dtype('complex128')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 框架介绍"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QPP 主要依赖一个名为**相位处理器**的量子线路来处理数据。对于偶数 $L \\in \\mathbb{N}$，我们定义一个 $n$-比特酉算子 $U$ 的相位处理器为\n",
    "\n",
    "$$\n",
    "V^L(U) := R_z^{(0)} R_y^{(0)} R_z^{(0)}\n",
    "\\left[ \\prod_{l=1}^{L/2}\n",
    "    \\begin{bmatrix}\n",
    "        U^\\dagger & 0 \\\\\n",
    "        0 & I^{\\otimes n}\n",
    "    \\end{bmatrix} R_y^{(0)} R_z^{(0)}\n",
    "    \\begin{bmatrix}\n",
    "        I^{\\otimes n} & 0 \\\\\n",
    "        0 & U\n",
    "    \\end{bmatrix} R_y^{(0)} R_z^{(0)}\n",
    "\\right]，\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "这里 $R_y^{(0)}$ 和 $R_z^{(0)}$ 为作用在处理器第一个比特上的旋转门，且其旋转角度取决于模拟的三角多项式。其具体的量子电路如下所示：\n",
    "\n",
    "![qpp-circuit](figures/QPP-fig-circuit.png \"图 1：相位处理器的量子实现（层数 L 为偶数）\")\n",
    "\n",
    "注：当层数 $L$ 为奇数时，我们定义相应的相位处理器为\n",
    "\n",
    "$$ \n",
    "V^{L}(U) = V^{L - 1}(U) \n",
    "    \\begin{bmatrix}\n",
    "        U^\\dagger & 0 \\\\\n",
    "        0 & I^{\\otimes n}\n",
    "    \\end{bmatrix} R_y^{(0)} R_z^{(0)} \\tag{2}。\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QPP 有两个核心功能：**相位演化**与**相位提取**。 具体来说，假设一个有如下谱分解的 $n$-（量子）比特酉矩阵 \n",
    "\n",
    "$$\n",
    "U = \\sum_{j=0}^{2^n - 1} e^{i \\tau_j} |\\chi_j \\rangle \\langle \\chi_j |。 \\tag{3}\n",
    "$$\n",
    "\n",
    "那么论文 [[1]](https://arxiv.org/abs/2209.14278) 中的 **Theorem 5** 证明，对于任意满足 $||\\textbf{c}||_1 \\leq 1$ 的（复数）三角多项式 $F(x) = \\sum_{j = -L}^L c_j e^{ijx}$，和任意量子态 $| \\psi \\rangle = \\sum_{j=0}^{2^n - 1} \\alpha_j |\\chi_j \\rangle$，我们都可以找到一个层数为 $2L$ 的相位处理器 $V^{2L}(U)$ 使得\n",
    "\n",
    "$$\n",
    "\\left( \\langle 0 | \\otimes I^{\\otimes n} \\right) V^{2L}(U) |0, \\psi \\rangle = \\sum_{j=0}^{2^n - 1} \\alpha_j F(\\tau_j) | \\chi_j \\rangle。 \\tag{4}\n",
    "$$\n",
    "\n",
    "进一步地，若 $F$ 的值域为实数域，那么对于任意一个量子态 $\\rho$, 论文 [[1]](https://arxiv.org/abs/2209.14278) 中的 **Theorem 6** 证明，我们可以找到一个层数为 $L$ 的相位处理器 $V^{L}(U)$ 使得\n",
    "\n",
    "$$\n",
    "\\text{Tr} \\left[ Z^{(0)} \\cdot V^{L}(U) \\rho V^{L}(U)^\\dagger \\right] = \\sum_{j = 0}^{2^n - 1} p_j F(\\tau_j)， \\tag{5}\n",
    "$$\n",
    "\n",
    "这里 $p_j = \\langle \\chi_j | \\rho | \\chi_j \\rangle$，$Z^{(0)}$ 是作用在第一个比特上的泡利-$Z$ 可观测量。下面我们将用这两个功能来完成酉矩阵的相位搜索、哈密顿量的实时演化和量子态的函数变换的迹估计。\n",
    "\n",
    "注：由于哈密顿量和量子态不是酉矩阵，无法被 QPP 直接处理，因此我们需要使用**比特化块编码**（qubitized block encoding，下文统一称作块编码）来拿到非酉矩阵的本征信息。\n",
    "另外，比特化块编码的本征相位和其编码数据的本征值还存在一个 $\\arccos$ 的关系。所以对于一个函数 $f$ 的变换，QPP 实际需要模拟的函数是 $F(x) = f(\\cos(x))$。\n",
    "更多块编码的介绍详见[量子信号处理与量子奇异值变换](https://qml.baidu.com/tutorials/quantum-simulation/quantum-signal-processing-and-quantum-singular-value-transformation.html)教程；比特化的知识可以参考论文 [[1]](https://arxiv.org/abs/2209.14278) 和论文 [[3]](https://quantum-journal.org/papers/q-2019-07-12-163/)。\n",
    "\n",
    "以下是本教程的环境设置。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_qubits = 3 # 各问题的主寄存器大小\n",
    "num_block_qubits = num_qubits + 1 # 使用 block encoding 的辅助比特数\n",
    "aux_qubits = list(range(num_block_qubits + 1)) # 辅助寄存器所在的比特索引\n",
    "sys_qubits = list(range(num_block_qubits + 1, num_block_qubits + 1 + num_qubits)) # 主寄存器所在的比特索引"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 应用：量子相位搜索"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "求解一个酉算子的本征相位是量子计算的核心问题之一，这个问题也叫做量子相位估计（quantum phase estimation，简称 QPE）。该问题的设定为：已知一个酉算子 $U$ 和它的一个本征态 $| \\psi \\rangle$，求解该本征态所对应的本征相位；进一步地，若 $| \\psi \\rangle$ 不再是 $U$ 的本征态，则转而求解（任意）一个与 $| \\psi \\rangle$ 相交的本征态所对应的本征相位。\n",
    "\n",
    "以下是该问题的实验设定。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "pq.set_backend('state_vector') # 切换态矢量后端\n",
    "\n",
    "U = unitary_random(num_qubits) # 酉算子\n",
    "psi = random_state(num_qubits) # 输入态（矢量）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QPP 可以模拟阶梯函数即\n",
    "\n",
    "$$\n",
    "f(x) = \\begin{cases}\n",
    "    1 & \\textrm{if }\\, x \\geq 0 \\\\\n",
    "    0 & \\textrm{if }\\, x < 0\n",
    "\\end{cases}。 \\tag{6}\n",
    "$$\n",
    "\n",
    "QPP 通过测量辅助比特的方式来二分搜索到一个本征相位所在的小区间；然后将这个区间放大再重复操作，从而完成本征相位的估计任务。我们称该算法为**量子相位搜索**（quantum phase search，简称 QPS）算法。相较于传统 QPE 算法，QPS 算法可以在同等精度和资源下获得成功概率的指数提升。QPS 算法的具体细节可以移步至论文 [[1]](https://arxiv.org/abs/2209.14278) 的第二节。\n",
    "\n",
    "量桨的 QPP 模块有一个内置函数 `qps`，其实现了完整的 QPS 算法。运行结束后，我们可以将结果与理论值对比来验证算法的正确性。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computations of angles for QPP are completed with mean error 1.692667280042757e-07\n",
      "估计特征相位与输出态特征相位的匹配误差为 2.104050267612231e-11\n",
      "输出态与输入态的 overlap 为 0.12645570235623324\n"
     ]
    }
   ],
   "source": [
    "# 使用 QPS 算法获取一个本征相位和主系统的输出态\n",
    "phase_estimate, output_state = qps(U, psi)\n",
    "\n",
    "# 拿到输出态对应的本征相位\n",
    "phase_expect = np.log((output_state.bra @ U @ output_state.ket).item()) / 1j\n",
    "\n",
    "print(f\"估计本征相位与输出态本征相位的匹配误差为 {np.abs(phase_expect - phase_estimate)}\")\n",
    "print(f\"输出态与输入态的 overlap 为 {abs_norm(output_state.bra @ psi.ket)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上可见，虽然输入态 $| \\psi \\rangle$ 不是 $U$ 的本征态，我们仍可以拿到与 $| \\psi \\rangle$ 内积不为 $0$ 的本征态对应的本征相位。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 应用：哈密顿量模拟"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一个随着时间演化的 $n$-比特量子系统可以由一个哈密顿量 $H$ 和这个量子系统的初始态 $\\rho$ 决定。该系统在时间 $t$ 的量子态可以表示为 $\\rho_t = e^{-iHt} \\rho e^{iHt}$，这里 $e^{-iHt}$ 称为该系统在时刻 $t$ 下的演化算子。那么哈密顿量模拟的问题设定为：已知一个量子系统的哈密顿量的块编码 $U$ 和初始态 $\\rho$，近似制备该系统在时刻 $t$ 下的量子态 $\\rho_t$。\n",
    "\n",
    "以下是该问题的实验设定。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pq.set_backend(\"density_matrix\") # 切换密度矩阵后端\n",
    "\n",
    "H = hermitian_random(num_qubits) # 量子系统的哈密顿量\n",
    "U = block_enc_herm(H, num_block_qubits) # 哈密顿量的比特化块编码\n",
    "rho = random_state(num_qubits) # 初始态（密度矩阵）\n",
    "t = 9 # 演化时间\n",
    "L = 40 # 相位处理器的层数，即模拟演化函数的精度"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "量桨 QPP 模块的 `hamiltonian_laurent` 可以通过 [Jacobi-Anger 展开](https://en.wikipedia.org/wiki/Jacobi%E2%80%93Anger_expansion) 提供演化函数的近似三角多项式 $P$；然后用 `Q_generation` 算出该多项式的余式 $Q$（两者满足 $PP^* + QQ^* = 1$）；函数 `qpp_angle_approximator` 负责估算相位处理器的所有旋转门角度；最后函数 `qpp_cir` 完成相位处理器的构建。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computations of angles for QPP are completed with mean error 1.1413227206232605e-07\n"
     ]
    }
   ],
   "source": [
    "# 准备近似模拟演化函数的三角多项式和其余式, 这里乘个略小于 1 的数来保证余式的成功计算\n",
    "P = hamiltonian_laurent(-t, L) * 0.999999\n",
    "Q = Q_generation(P)\n",
    "\n",
    "# 计算旋转角度，其中 theta 对应 Ry 门，phi 对应 Rz 门\n",
    "list_theta, list_phi = qpp_angle_approximator(P, Q)\n",
    "\n",
    "cir = qpp_cir(list_theta, list_phi, U) # 构建相位处理器\n",
    "cir.collapse(aux_qubits, desired_result=0, if_print=True) # 在线路末端添加解码（测量）操作\n",
    "input_state = zero_state(num_block_qubits + 1).kron(rho) # 准备输入态，其中辅助比特的输入态为 |0><0|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "qubits [0, 1, 2, 3, 4] collapse to the state |00000> with probability 0.9999974132814917\n",
      "输出态与期望态的迹距离为 4.903932274526325e-07\n"
     ]
    }
   ],
   "source": [
    "# 拿到输出态并移除辅助比特，将输出态与理论值对比\n",
    "output_state = partial_trace_discontiguous(cir(input_state), preserve_qubits=sys_qubits)\n",
    "rho.evolve(H, t)\n",
    "print(f\"输出态与期望态的迹距离为 {trace_distance(output_state, rho).item()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上可见，QPP 通过相位演化，将一个 $H$ 的块编码转换为 $e^{-iHt}$ 的块编码，从而完成量子态 $\\rho_t$ 的成功制备。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 应用：量子态熵估计"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如何计算量子态函数变换后的迹是量子熵估计的核心问题。该问题的数学设定为：已知两个量子态 $\\rho$、 $\\sigma$ 和一个函数 $f:~\\mathbb{R}_+~\\to~\\mathbb{R}$，估计 $\\text{Tr}\\left[ \\rho f \\left( \\sigma \\right) \\right]$。\n",
    "\n",
    "注：量子态 $\\sigma$ 的 $f$ 变换定义为 \n",
    "\n",
    "$$ \n",
    "f(\\sigma) = \\sum_j f(q_j) | \\psi_j \\rangle \\langle \\psi_j |， \\tag{7}\n",
    "$$\n",
    "这里 $\\{q_j\\}$ 和 $\\{ | \\psi_j \\rangle \\}$ 分别为 $\\sigma$ 的本征值和本征态。 \n",
    "\n",
    "可以看到，若 $f$ 为幂函数，该问题可以用来解决量子瑞丽熵和瑞丽散度的估计；而当 $f(x) = \\log(x)$ 时，该问题则转化为量子冯诺依曼熵或相对熵的估计问题。\n",
    "\n",
    "QPP 的相位提取功能可以有效地解决迹估计问题，从而进一步地完成各类量子熵的近似计算。用于处理量子态的 QPP 线路如图所示：\n",
    "\n",
    "![qpp](figures/QPP-fig-state.png \"图 2：用于迹估计的相位处理器（层数为偶数），这里 m 为用于块编码的辅助比特数。\")\n",
    "\n",
    "这里作用在 AB 系统上的 $U_\\rho$ 为 $\\rho$ 的**纯化模型**，满足\n",
    "\n",
    "$$\n",
    "\\text{Tr}_B \\left[ U_\\rho \\left( | 0 \\rangle_A \\langle 0 |_A \\otimes | 0 \\rangle_B \\langle 0 |_B \\right) U_\\rho^\\dagger \\right] = \\rho。 \\tag{8}\n",
    "$$\n",
    "\n",
    "该模型可以用量子电路实现，所以在量子熵的研究中被广泛使用。$U_\\sigma$ 也为 $\\sigma$ 的纯化模型。不同的是，为了拿到 $\\sigma$ 的本征信息，QPP 需要将 $U_\\sigma$ 进一步转为 $\\sigma$ 的块编码 $\\widehat{U}_\\sigma$。更多内容详见论文 [[1]](https://arxiv.org/abs/2209.14278) 的第四节。\n",
    "\n",
    "我们将使用 QPP 模块来估计 $\\text{Tr} \\left[ \\rho \\sigma^{\\alpha - 1} \\right]$。以下是该问题的实验设定。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pq.set_backend(\"density_matrix\") # 切换密度矩阵后端\n",
    "\n",
    "rho = random_state(num_qubits) # 输入态（密度矩阵），这里我们不关注上图的 B 系统\n",
    "U_sigma_hat = purification_block_enc(num_qubits, num_block_qubits) # 构建一个纯化模型的比特化块编码\n",
    "sigma = State(U_sigma_hat[:2**num_qubits, :2**num_qubits]) # 通过块编码获取随机 sigma\n",
    "assert is_density_matrix(sigma) == (True, num_qubits)\n",
    "\n",
    "alpha = np.random.rand() * 4 + 1 # 在 [1, 5) 之间随机选择 alpha\n",
    "H = Hamiltonian([(1.0, \"z0\")]) # 第一个比特上的泡利-Z 可观测量\n",
    "input_state = zero_state(num_block_qubits + 1).kron(rho) # 相位处理器的输入态"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "量桨 QPP 模块里的 `simulation_cir` 函数可以为 $f$ 设计适应性的相位处理器。具体地，通过找到函数 $f$ 的近似三角多项式 $F$，QPP 模块可以在机械误差下将 $F$ 转为相位处理器需要模拟的三角多项式，然后通过估算第一个比特上泡利-$Z$ 可观测量的期望值来获取最终的结果。我们可以借助 `paddle_quantum.linalg.herm_transform` 来评估模拟的精准度。\n",
    "\n",
    "注：函数 $f$ 需要满足 $f\\left(\\left[0, 1\\right]\\right) \\subseteq \\left[-1, 1\\right]$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computations of angles for QPP are completed with mean error 2.073248366556774e-07\n",
      "模拟 Tr[rho * sigma^3.522308957539911] 的误差为 2.265239344059755e-05\n"
     ]
    }
   ],
   "source": [
    "# 估计 Tr[rho * sigma^(alpha - 1)]\n",
    "cir = simulation_cir(lambda x: (np.cos(x) ** 2) ** ((alpha - 1) / 2), U_sigma_hat)\n",
    "val = cir(input_state).expec_val(H).item()\n",
    "expect_val = paddle.trace(rho.data @ herm_transform(lambda x: x ** (alpha - 1), sigma)).real().item()\n",
    "print(f\"模拟 Tr[rho * sigma^{alpha - 1}] 的误差为 {np.abs(val - expect_val)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上可见，QPP 通过相位提取来完成 $\\text{Tr} \\left[ \\rho \\sigma^{\\alpha - 1} \\right]$ 的精确估计。因此，当 $\\rho = \\sigma$ 时，QPP 就可以估算 $\\rho$ 的 $\\alpha$-阶量子瑞丽熵 \n",
    "\n",
    "$$\n",
    "S_\\alpha(\\rho) = \\frac{1}{1 - \\alpha}\\log \\text{Tr} \\left( \\rho^{\\alpha } \\right)。  \\tag{9}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 小结\n",
    "\n",
    "通过上述三个应用，本教程展示了 QPP 框架在酉算子、哈密顿量和量子态相关问题上的强大处理能力，以及如何用量桨的 QPP 模块去完成框架的构建和运算。我们也期望使用 QPP 框架去解决更多问题包括但不限于量子蒙特卡洛问题、酉矩阵迹估计问题以及机器学习问题。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "## 参考文献\n",
    "\n",
    "[1] Wang, Xin, et al. \"Quantum Phase Processing: Transform and Extract Eigen-Information of Quantum Systems.\" [arXiv preprint arXiv:2209.14278 (2022).](https://arxiv.org/abs/2209.14278)\n",
    "\n",
    "[2] Yu, Zhan, et al. \"Power and limitations of single-qubit native quantum neural networks.\" [arXiv preprint arXiv:2205.07848 (2022).](https://arxiv.org/abs/2205.07848)\n",
    "\n",
    "[3] Low, Guang Hao, and Isaac L. Chuang. \"Hamiltonian simulation by qubitization.\" [Quantum 3 (2019): 163.](https://quantum-journal.org/papers/q-2019-07-12-163/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.13 ('pq')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "08942b1340a5932ff3a93f52933a99b0e263568f3aace1d262ffa4d9a0f2da31"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
